parameters

if(is.null(input$cluster)) {
  input$cluster = '0'
  print("Warning: setting cluster to 0")
}
if(is.null(input$cluster5)) {
  input$cluster5 = '0'
  print("Warning: setting cluster5 to 0")
}
if(is.null(input$clusters)) {
  input$clusters = '0'
  print("Warning: setting clusters to 0")
}

## [1] "Warning: setting clusters to 0"

if(is.null(input$clusters1)) {
  input$clusters1 = '0'
  print("Warning: setting clusters1 to All")
}
if(is.null(input$clusters2)) {
  input$clusters2 = '0'
  print("Warning: setting clusters2 to 0")
}
if(is.null(input$clusters3)) {
  input$clusters3 = '0'
  print("Warning: setting clusters4 to 0")
}
if(is.null(input$clusters4)) {
  input$clusters4 = 'All'
  print("Warning: setting clusters4 to All")
}

dataTables = list()
dataTables$log2cpmOrg <- log2cpm
dataTables$tsne.data <- tsne.data
dataTables$tsne.dataOrg <- tsne.data
featuredata <-
  featuredata[which(featuredata$Chromosome.Name %in% c(unlist(lapply(
    seq(1, 22, 1), toString
  )), c("X", "Y", "MT", "N"))),]
featuredata$Associated.Gene.Name <-
  toupper(featuredata$Associated.Gene.Name)
featuredata <- featuredata[rownames(log2cpm), ]
dataTables$featuredataOrg <- featuredata
dataTables$positiveCells <- NULL
dataTables$positiveCellsAll <- NULL

# take only genes that are in all tables
rnames = rownames(featuredata)
rnames = rnames[rnames %in% rownames(log2cpm)]
rnames = rnames[rnames %in% rownames(gbm)]

dataTables$log2cpm <- log2cpm[rnames, ]
dataTables$gbm = gbm[rnames, ]
dataTables$gbm_log = gbm_log[rnames, ]
dataTables$featuredata <- featuredata[rnames, ]

geneC = colSums(log2cpm>0, na.rm = TRUE)
medianENSG = median(t(geneC))

umiC = colSums(log2cpm, na.rm = TRUE)
medianUMI = median(t(umiC))

minG = input$minGenes
gbm = as.matrix(exprs(dataTables$gbm))
goodCols = colSums(gbm) > minG

maxG = input$maxGenes
goodCols = goodCols & ( colSums(gbm) <= maxG)

geneNames = input$minExpGenes
ids = which(grepl(geneNames, dataTables$featuredata$Associated.Gene.Name))
goodCols = goodCols & (colSums(gbm[ids, ]) > 0)


useCells = goodCols

#genes to be explicitly removed
ipIDs = input$selectIds
if(nchar(ipIDs)>0) {
  rmIds = !grepl(input$selectIds, dataTables$featuredata$Associated.Gene.Name)
}else 
{
  rmIds = rep(TRUE,nrow(dataTables$log2cpm))
}

# gene groups to be included
if(!is.null(input$geneListSelection)){
  selectedgeneList = get_selected(input$geneListSelection)
  if(length(selectedgeneList)>0){
    selGenes = c()
    for(sIdx in 1:length(selectedgeneList)){
      print(sIdx)
      att = attr(selectedgeneList[[sIdx]], "ancestry")
      if(length(att)>0){
        selGenes = c(selGenes , geneLists[[att]][[selectedgeneList[[sIdx]]]])
      }
    }
    selGenes = unique(selGenes)
    rmIds = rownames(dataTables$log2cpm) %in% selGenes & rmIds
  }
}

# overall gene expression Min
minGene <- input$minGenesGS
if(!is.null(minGene)){
  selGenes = rowSums(as.matrix(exprs(dataTables$gbm[,useCells]))) >=minGene
  rmIds = rmIds & selGenes
}

if(DEBUG)cat(file=stderr(), "useGenes: done\n")
useGenes = rmIds

featureDataReact = dataTables$featuredata[useGenes, ]
featureData = dataTables$featuredata[useGenes, ]

if(DEBUG)cat(file=stderr(), "gbm\n")
gbm = dataTables$gbm[useGenes, useCells]
if(DEBUG)cat(file=stderr(), "gbm:DONE\n")

if(DEBUG)cat(file=stderr(), "gbm_log\n")
gbm_log = dataTables$gbm_log[useGenes, useCells]
if(DEBUG)cat(file=stderr(), "gbm_log:Done\n")

if(DEBUG)cat(file=stderr(), "scaterReads\n")
fd = featureDataReact

counts = as.matrix(exprs( gbm))

anno = pData(gbm)
anno$sample_id = anno$barcode
anno$fixed = "red"
# anno$individual= "NA1"
# anno$replicate = "r1"
# anno$well = "A01"
# anno$batch = "b1"
pheno_data <- new("AnnotatedDataFrame", anno)
rownames(pheno_data) <- pheno_data$sample_id

reads = as.matrix(counts)
rownames(reads) = make.unique(fd[rownames(reads),"Associated.Gene.Name"])
rownames(reads)[is.na(rownames(reads)) ] = "na"
reads = SingleCellExperiment(
  assays = list(counts = reads), 
  colData = anno
)
ercc <- rownames(reads)[grepl("ERCC-", rownames(reads))]

mt = rownames(reads)[grepl("^MT",rownames(reads))]

reads <- scater::calculateQCMetrics(
  reads
)
filter_by_expr_features <- (reads$total_features > 200)
reads$use <- (
  # sufficient features (genes)
  filter_by_expr_features 
  # sufficient molecules counted
  #filter_by_total_counts &
  # sufficient endogenous RNA
  #filter_by_ERCC &
  # remove cells with unusual number of reads in MT genes
  #filter_by_MT
)


scaterReads = reads

if(DEBUG)cat(file=stderr(), "log2cpm\n")
if(DEBUG)cat(file=stderr(), "log2cpm:Done\n")

log2cpm = dataTables$log2cpm[useGenes, useCells]

if(DEBUG)cat(file=stderr(), "pca\n")

pca = run_pca(gbm_log)

## Variance explained by PCs: 0.02963652

if(DEBUG)cat(file=stderr(), "tsne\n")

set.seed(seed = seed)
tsne = run_tsne(pca, dims = 3, perplexity = 30, theta = 0.5)

## Using number of of pcs: 10

if(DEBUG)cat(file=stderr(), "tsne: done\n")

if(DEBUG)cat(file=stderr(), "kmClustering\n")
clustering=list()

kNr = 10
for(kNr in 2:10) {
  set.seed(seed = seed)
  km = run_kmeans_clustering(pca, k=kNr)
  clustering[[paste0("kmeans_",kNr,"_clusters")]] = data.frame("Barcode" = rownames(data.frame(km$cluster)), "Cluster" = km$cluster)
}

## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10 
## Using number of of pcs: 10

kmClustering = clustering
if(DEBUG)cat(file=stderr(), "kmClustering:done\n")

if(DEBUG)cat(file=stderr(), "tsne.data\n")

tsne.data = data.frame(tsne$Y)
colnames(tsne.data) = c("V1", "V2", "V3")
tsne.data$dbCluster = clustering$kmeans_10_clusters$Cluster-1
rownames(tsne.data) = clustering$kmeans_10_clusters$Barcode
if(DEBUG)cat(file=stderr(), "tsne.data: done\n")

set.seed(seed = seed)
prioritized_genes =     prioritize_top_genes(gbm, tsne.data$dbCluster, "sseq",
                                             logscale = FALSE, 
                                             min_mean=0.5, 
                                             p_cutoff=0.05,
                                             order_by='pvalue')

## Computing differential expression parameters...
## Comparing Cluster 0 against the other clusters...
## Comparing Cluster 1 against the other clusters...
## Comparing Cluster 2 against the other clusters...
## Comparing Cluster 3 against the other clusters...
## Comparing Cluster 4 against the other clusters...
## Comparing Cluster 5 against the other clusters...
## Comparing Cluster 6 against the other clusters...
## Comparing Cluster 7 against the other clusters...
## Comparing Cluster 8 against the other clusters...
## Comparing Cluster 9 against the other clusters...

if(DEBUG)cat(file=stderr(), "dge\n")
if(!is.null(input$db1) & ! is.null(input$db2)){
  subsetData <- subset(tsne.data, dbCluster == input$clusters1)
  cells.1 <- rownames(brushedPoints(subsetData, input$db1))
  cells.2 <- rownames(brushedPoints(subsetData, input$db2))

  subsetExpression <- log2cpm[, union(cells.1, cells.2)]
  
  genes.use <- rownames(subsetExpression)
  data.1 = apply(subsetExpression[genes.use, cells.1], 1, expMean)
  data.2 = apply(subsetExpression[genes.use, cells.2], 1, expMean)
  total.diff = (data.1 - data.2)
  
  genes.diff = names(which(abs(total.diff) > .2))
  genes.use = ainb(genes.use, genes.diff)
  
  toReturn <-
    DiffExpTest(subsetExpression, cells.1, cells.2, genes.use = genes.use)
  toReturn[, "avg_diff"] = total.diff[rownames(toReturn)]
  toReturn$Associated.Gene.Name <-
    featureData[rownames(toReturn), 'Associated.Gene.Name']
  selectedDge <- toReturn
  cat(stderr(), rownames(toReturn)[1:5])
  
  dge <-     toReturn
  if(DEBUG)cat(file=stderr(), "dge:done11\n")
}else{
  dge=NA
}

## 2 ENSG00000186891 ENSG00000116209 ENSG00000031698 ENSG00000265241 ENSG00000163191

# selectedDge <- reactiveValues()

p <-
  plot_ly(
    tsne.data,
    x = ~ V1,
    y = ~ V2,
    z = ~ V3,
    type = "scatter3d",
    color =  ~ dbCluster,
    hoverinfo = "text",
    text = paste('Cluster:', tsne.data$dbCluster),
    mode = 'markers',
    marker =
      list(
        line = list(width = 0),
        size = rep(10, nrow(tsne.data)),
        sizeref = 3
      )
  )
if(DEBUG)cat(file=stderr(), "output$tsne_main: done\n")
layout(p)

line1<-paste('No. of cells:', dim(log2cpm)[2],sep='\t')
line2<-paste('Median UMIs:', medianUMI,sep='\t')
line3<-paste('Median Genes:', medianENSG,sep='\t')
line5<-paste('No. of reads:',  dim(log2cpm)[1],sep='\t')
HTML(
  paste0("Summary statistics of this dataset:", '<br/>','<br/>',
         line1, '<br/>', 
         line2, '<br/>',
         line3, '<br/>',
         line5
  )
)

Summary statistics of this dataset:

No. of cells: 904
Median UMIs: 2231.12313456375
Median Genes: 1216
No. of reads: 14598

hist(colSums(as.matrix(exprs(gbm))), breaks = 50, main="histogram of number of UMI per cell")

barplot(pca$var_pcs, main="Variance captured by first PCs")

barplot(pca$var_pcs, main="Variance captured by first PCs")

gbmMat = as.matrix(exprs(dataTables$gbm))
fd = dataTables$featuredata
dt = fd[useGenes,c("Associated.Gene.Name", "Gene.Biotype", "Description")]
dt$rowSums = rowSums(gbmMat[useGenes,useCells])
dt$rowSamples = rowSums(gbmMat[useGenes,useCells]>0)
DT::datatable(dt)

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

gbmMat = as.matrix(exprs(dataTables$gbm))
fd = dataTables$featuredata
dt = fd[useGenes,c("Associated.Gene.Name", "Gene.Biotype", "Description")]
dt$rowSums = rowSums(gbmMat[useGenes,useCells])
dt$rowSamples = rowSums(gbmMat[useGenes,useCells]>0)
DT::datatable(dt)

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

if(DEBUG)cat(file=stderr(), "output$tsne_plt\n")
geneid <- rownames(featureData[which(featureData$Associated.Gene.Name ==
                                       toupper(input$gene_id)), ])[1]

expression <- log2cpm[geneid, ]

tsneData <- cbind(tsne.data, t(expression))
names(tsneData)[names(tsneData) == geneid] <- 'values'

p <-
  plot_ly(
    tsneData,
    x = ~ V1,
    y = ~ V2,
    z = ~ V3,
    type = "scatter3d",
    hoverinfo = "text",
    text = paste('Cluster:', tsne.data$dbCluster),
    mode = 'markers',
    marker = list(
      size = 2,
      line = list(width = 0),
      color =  ~ values,
      colors = 'Greens'
    )
  )
layout(p, title = toupper(input$gene_id))

if(DEBUG)cat(file=stderr(), "output$clusterPlot\n")
geneid <- rownames(featureData[which(featureData$Associated.Gene.Name ==
                                       toupper(input$gene_id)), ])[1]

expression <- log2cpm[geneid, ]


tsneData <- cbind(tsne.data, t(expression))
names(tsneData)[names(tsneData) == geneid] <- 'values'


if(is.null(input$cluster)) input$cluster = 0
subsetData <- subset(tsneData, dbCluster == input$cluster)

## Warning in dbCluster == input$cluster: longer object length is not a
## multiple of shorter object length

p1 <-
  ggplot(subsetData,
         aes_string(x = input$dimension_x, y = input$dimension_y)) +
  geom_point(aes_string(size = 2, color = 'values')) +
  geom_point(shape = 1,
             size = 4,
             colour = "black") +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 90,
      size = 12,
      vjust = 0.5
    ),
    axis.text.y = element_text(size = 12),
    strip.text.x = element_text(size = 16),
    strip.text.y = element_text(size = 14),
    axis.title.x = element_text(face = "bold", size = 16),
    axis.title.y = element_text(face = "bold", size = 16),
    legend.position = "none"
  ) +
  ggtitle(paste(toupper(input$gene_id), input$cluster, sep = '-Cluster')) +
  scale_colour_gradient2(low = 'grey50', high = "red")
p1

if(DEBUG)cat(file=stderr(), "output$gene_vio_plot\n")

geneid <- rownames(featureData[which(featureData$Associated.Gene.Name ==
                                       toupper(input$gene_id)), ])[1]

expression <- log2cpm[geneid, ]

validate(need(is.na(sum(expression)) != TRUE, ''))

tsneData <- cbind(tsne.data, t(expression))
names(tsneData)[names(tsneData) == geneid] <- 'values'
#tsne.data<-subset(tsne.data,dbCluster!=0)

p1 <-
  ggplot(tsneData, aes(factor(dbCluster), values, fill = factor(dbCluster))) +
  geom_violin(scale = "width") +
  stat_summary(
    fun.y = median,
    geom = "point",
    size = 5,
    color = 'black'
  ) +
  stat_summary(fun.data = n_fun, geom = "text") +
  theme_bw() +
  theme(
    axis.text.x = element_text(
      angle = 90,
      size = 12,
      vjust = 0.5
    ),
    axis.text.y = element_text(size = 12),
    strip.text.x = element_text(size = 16),
    strip.text.y = element_text(size = 14),
    axis.title.x = element_text(face = "bold", size = 16),
    axis.title.y = element_text(face = "bold", size = 16),
    legend.position = "none"
  ) +
  xlab('Cluster') +
  ylab('Expression') +
  ggtitle(toupper(input$gene_id))
if(DEBUG)cat(file=stderr(), "output$gene_vio_plot:done\n")
p1

if(DEBUG)cat(file=stderr(), "output$panelPlot\n")

genesin <- input$panelplotids
genesin <- toupper(genesin)
genesin <- gsub(" ", "", genesin, fixed = TRUE)
genesin <- strsplit(genesin, ',')
genesin<-genesin[[1]]

par(mfrow=c(ceiling(length(genesin)/4),4), mai = c(0, 0., 0., 0.))
rbPal <- colorRampPalette(c('#f0f0f0','red'))

if (input$clusters4 == 'All') 
{
  for (i in 1:length(genesin)){
    Col <- rbPal(10)[
      as.numeric(
        cut(
          as.numeric(
            log2cpm[
              rownames(featureData[which(featureData$Associated.Gene.Name==genesin[i]),])
              ,]
          ),breaks = 10))]
    plot(tsne.data[,input$dimension_x4],tsne.data[,input$dimension_y4],col=Col,pch=16,axes = FALSE,frame.plot = TRUE, ann=FALSE)
    title(genesin[i],line=-1.2,adj = 0.05,cex.main=2)
  }
}else{
  for (i in 1:length(genesin)){
    
    subsetTSNE <- subset(tsne.data, dbCluster == input$clusters4)
    
    Col <- rbPal(10)[
      as.numeric(
        cut(
          as.numeric(
            log2cpm[
              rownames(featureData[which(featureData$Associated.Gene.Name==genesin[i]),])
              ,]
          ),breaks = 10))]
    
    names(Col)<-rownames(tsne.data)
    plotCol<-Col[rownames(subsetTSNE)]
    plot(subsetTSNE[,input$dimension_x4],subsetTSNE[,input$dimension_y4],col=plotCol,pch=16,axes = FALSE,frame.plot = TRUE, ann=FALSE)
    title(genesin[i],line=-1.2,adj = 0.05,cex.main=2)
    
  }
}

scater::plotQC(scaterReads, type = "highest-expression", col_by_variable="fixed")

genesin <- input$heatmap_geneids
genesin <- toupper(genesin)
genesin <- gsub(" ", "", genesin, fixed = TRUE)
genesin <- strsplit(genesin, ',')

map <- rownames(featureData[which(featureData$Associated.Gene.Name %in% genesin[[1]]), ])
cat(file = stderr(), length(map))

expression <- log2cpm[map, ]

if(is.na(sum(expression)) == TRUE){
  print('Gene symbol incorrect or genes not expressed')
}else{
  
  tsne.data <- tsne.data[order(tsne.data$dbCluster), ]
  
  expression <- expression[, rownames(tsne.data)]
  expression <- expression[complete.cases(expression), ]
  
  annotation <- data.frame(factor(tsne.data$dbCluster))
  rownames(annotation) <- colnames(expression)
  colnames(annotation) <- c('Cluster')
  
  h <-
    pheatmap(
      as.matrix(expression),
      cluster_rows = TRUE,
      cluster_cols = FALSE,
      scale = 'row',
      fontsize_row = 10,
      labels_col = colnames(expression),
      labels_row = featureData[rownames(expression), 'Associated.Gene.Name'],
      show_rownames = TRUE,
      annotation_col = annotation,
      show_colnames = FALSE,
      annotation_legend = TRUE,
      breaks = seq(-6, 6, by = .12),
      colorRampPalette(rev(brewer.pal(
        n = 6, name =
          "RdBu"
      )))(100)
      
    )
  print(h)
}

## $tree_row
## 
## Call:
## hclust(d = d, method = method)
## 
## Cluster method   : complete 
## Distance         : euclidean 
## Number of objects: 4 
## 
## 
## $tree_col
## [1] NA
## 
## $kmeans
## [1] NA
## 
## $gtable
## TableGrob (5 x 6) "layout": 7 grobs
##   z     cells                 name                          grob
## 1 1 (4-4,1-1)             row_tree polyline[GRID.polyline.16057]
## 2 2 (4-4,3-3)               matrix       gTree[GRID.gTree.16059]
## 3 3 (4-4,4-4)            row_names         text[GRID.text.16060]
## 4 4 (3-3,3-3)       col_annotation         rect[GRID.rect.16061]
## 5 5 (3-3,4-4) col_annotation_names         text[GRID.text.16062]
## 6 6 (3-5,6-6)    annotation_legend       gTree[GRID.gTree.16067]
## 7 7 (3-5,5-5)               legend       gTree[GRID.gTree.16070]

if(DEBUG)cat(file=stderr(), "output$clusterPlot2\n")
geneid <- rownames(featureData[which(featureData$Associated.Gene.Name ==
                                       toupper(input$gene_id_sch)), ])[1]

expression <- log2cpm[geneid, ]

validate(need(is.na(sum(expression)) != TRUE, ''))

tsneData <- cbind(tsne.data, t(expression))
names(tsneData)[names(tsneData) == geneid] <- 'values'

# if(DEBUG)cat(file=stderr(), paste("output$dge_plot1:---",input$clusters2,"---\n"))
if(!is.null(input$dimension_x2) & !is.null(input$dimension_y2)){
  subsetData <- subset(tsneData, dbCluster == input$clusters2)
  p1 <-
    ggplot(subsetData,
           aes_string(x = input$dimension_x2, y = input$dimension_y2)) +
    geom_point(aes_string(size = 2, color = 'values')) +
    geom_point(shape = 1,
               size = 4,
               colour = "black") +
    theme_bw() +
    theme(
      axis.text.x = element_text(
        angle = 90,
        size = 12,
        vjust = 0.5
      ),
      axis.text.y = element_text(size = 10),
      strip.text.x = element_text(size = 16),
      strip.text.y = element_text(size = 14),
      axis.title.x = element_text(face = "bold", size = 16),
      axis.title.y = element_text(face = "bold", size = 16),
      legend.position = "none"
    ) +
    ggtitle(paste(toupper(input$gene_id_sch), input$clusters2, sep =
                    '-Cluster')) +
    scale_colour_gradient2(low = 'grey50', high = "red")
  print(p1)
}

if(DEBUG)cat(file=stderr(), "output$selectedHeatmap\n")

genesin <- input$heatmap_geneids2
genesin <- toupper(genesin)
genesin <- strsplit(genesin, ',')
if(!is.null(input$scb1)){
  subsetData <-
    subset(tsne.data, dbCluster == input$clusters2)
  cells.1 <- rownames(brushedPoints(df=subsetData, 
                                    brush=input$scb1,
                                    xvar = input$scb1$mapping$x,
                                    yvar = input$scb1$mapping$y,
                                    panelvar1 = 
                                    ))
  map <- rownames(featureData[which(featureData$Associated.Gene.Name %in% genesin[[1]]), ])

  expression <- log2cpm[map, cells.1]

  expression <- expression[complete.cases(expression), ]
  cat(file = stderr(), rownames(expression))
  mColor <- max(expression)
  
  if( is.na(sum(expression)) == TRUE){
    print('Gene symbol incorrect or genes not expressed')
  }else{
    

    h <-
      pheatmap(
        as.matrix(expression),
        cluster_rows = TRUE,
        cluster_cols = TRUE,
        scale = 'row',
        fontsize_row = 10,
        labels_col = colnames(expression),
        labels_row = featureData[rownames(expression), 'Associated.Gene.Name'],
        show_rownames = TRUE,
        show_colnames = FALSE,
        breaks = seq(-6, 6, by = .12),
        colorRampPalette(rev(brewer.pal(
          n = 6, name =
            "RdBu"
        )))(100)
        
      )
    print(h)
  }
}

## $tree_row
## 
## Call:
## hclust(d = d, method = method)
## 
## Cluster method   : complete 
## Distance         : euclidean 
## Number of objects: 2 
## 
## 
## $tree_col
## 
## Call:
## hclust(d = d, method = method)
## 
## Cluster method   : complete 
## Distance         : euclidean 
## Number of objects: 44 
## 
## 
## $kmeans
## [1] NA
## 
## $gtable
## TableGrob (5 x 6) "layout": 5 grobs
##   z     cells      name                          grob
## 1 1 (2-2,3-3)  col_tree polyline[GRID.polyline.16121]
## 2 2 (4-4,1-1)  row_tree polyline[GRID.polyline.16122]
## 3 3 (4-4,3-3)    matrix       gTree[GRID.gTree.16124]
## 4 4 (4-4,4-4) row_names         text[GRID.text.16125]
## 5 5 (3-5,5-5)    legend       gTree[GRID.gTree.16128]

if(DEBUG)cat(file=stderr(), "output$plotCoExpression\n")

genesin <- input$mclustids
genesin <- toupper(genesin)
genesin <- strsplit(genesin, ',')

subsetData <-
  subset(tsne.data, dbCluster == input$clusters3)
cells.1 <- rownames(subsetData)


map <- rownames(featureData[which(featureData$Associated.Gene.Name %in% genesin[[1]]), ])

expression <- log2cpm[map, ]

if(is.na(sum(expression)) == TRUE){
  print ('Gene symbol incorrect or genes not expressed')
}else{
  
  bin <- expression
  bin[] <- 0
  
  for (i in 1:nrow(expression))
  {
    x <- Mclust(expression[i, ], G = 2)
    bin[i, ] <- x$classification
  }
  bin <- bin - 1
  allexprs <- apply(bin, 2, sum)
  plotexprs <- allexprs
  plotexprs[] <- 0
  plotexprs[allexprs >= length(rownames(bin))] <- 1
  positiveCells <- allexprs >= length(rownames(bin))
  positiveCellsAll <- plotexprs
    mergeExprs <- plotexprs[rownames(subsetData)]
  
  subsetData$CoExpression <- mergeExprs
  
  p1 <-
    ggplot(subsetData,
           aes_string(x = input$dimension_x3, y = input$dimension_y3)) +
    geom_point(aes_string(size = 2, color = 'CoExpression')) +
    geom_point(shape = 1,
               size = 4,
               colour = "black") +
    theme_bw() +
    theme(
      axis.text.x = element_text(
        angle = 90,
        size = 12,
        vjust = 0.5
      ),
      axis.text.y = element_text(size = 12),
      strip.text.x = element_text(size = 16),
      strip.text.y = element_text(size = 14),
      axis.title.x = element_text(face = "bold", size = 16),
      axis.title.y = element_text(face = "bold", size = 16),
      legend.position = "none"
    ) +
    #ggtitle(paste(toupper(input$gene_id),input$cluster,sep='-Cluster'))+
    scale_colour_gradient2(low = 'grey50', high = "red")
  print(p1)
}

## fitting ...
## 
  |                                                                       
  |                                                                 |   0%
  |                                                                       
  |======================                                           |  33%
  |                                                                       
  |===========================================                      |  67%
  |                                                                       
  |=================================================================| 100%

if(DEBUG)cat(file=stderr(), "output$onOffTable\n")

merge <- tsne.data
# if(DEBUG)cat(file=stderr(), paste("positiveCellsAll:---",positiveCellsAll,"---\n"))

merge$CoExpression <- positiveCellsAll
df <-
  as.data.frame(table(merge[, c('dbCluster', 'CoExpression')]))
dfOut <- cast(df, dbCluster ~ CoExpression)

## Using Freq as value column.  Use the value argument to cast to override this choice

colnames(dfOut) <- c("Cluster", 'OFF', 'ON')
rownames(dfOut) <- dfOut$Cluster
dfOut['Sum', ] <- c('', sum(dfOut$OFF), sum(dfOut$ON))

## Warning in `[<-.factor`(`*tmp*`, iseq, value = ""): invalid factor level,
## NA generated

DT::datatable(dfOut)

if(!is.null(input$dimension_x1) & !is.null(input$dimension_y1)){
  subsetData <- subset(tsne.data, dbCluster == input$clusters1)
  p1 <-
    ggplot(subsetData,
           aes_string(x = input$dimension_x1, y = input$dimension_y1)) +
    geom_point() +
    geom_point(shape = 1,
               size = 4,
               color = "black") +
    theme_bw() +
    theme(
      axis.text.x = element_text(
        angle = 90,
        size = 12,
        vjust = 0.5
      ),
      axis.text.y = element_text(size = 12),
      strip.text.x = element_text(size = 16),
      strip.text.y = element_text(size = 14),
      axis.title.x = element_text(face = "bold", size = 16),
      axis.title.y = element_text(face = "bold", size = 16),
      legend.position = "none"
    ) +
    ggtitle(input$clusters1)
  print(p1)
}

if(DEBUG)cat(file=stderr(), "output$dge_plot2\n")
if(!is.null(input$dimension_x1) & !is.null(input$dimension_y1)){
  subsetData <- subset(tsne.data, dbCluster == input$clusters1)
  p1 <-
    ggplot(subsetData,
           aes_string(x = input$dimension_x1, y = input$dimension_y1)) +
    geom_point() +
    geom_point(shape = 1,
               size = 4,
               color = "black") +
    theme_bw() +
    theme(
      axis.text.x = element_text(
        angle = 90,
        size = 12,
        vjust = 0.5
      ),
      axis.text.y = element_text(size = 12),
      strip.text.x = element_text(size = 16),
      strip.text.y = element_text(size = 14),
      axis.title.x = element_text(face = "bold", size = 16),
      axis.title.y = element_text(face = "bold", size = 16),
      legend.position = "none"
    ) +
    ggtitle(input$clusters1)
  print(p1)
}

if(DEBUG)cat(file=stderr(), "output$dge\n")

if(!is.na(dge)){
  top.genes <- dge
  top.genes$Associated.Gene.Name <-
    featureData[rownames(top.genes), 'Associated.Gene.Name']
  if (dim(top.genes)[1] > 1) {
    DT::datatable(top.genes,
                  options = list(
                    orderClasses = TRUE,
                    lengthMenu = c(10, 30, 50),
                    pageLength = 10
                  ))
  }
}

## Warning in if (!is.na(dge)) {: the condition has length > 1 and only the
## first element will be used

if(DEBUG)cat(file=stderr(), "output$crHeat_plot1\n")

example_K <- 10 
example_Cols <- rev(brewer.pal(10,"Set3")) # customize plotting colors

cells_to_plot <- order_cell_by_clusters(gbm, tsne.data$dbCluster)

example_col = example_Cols[1:example_K]
p = gbm_pheatmap(gbm=log_gene_bc_matrix(gbm), 
                 genes_to_plot=prioritized_genes, 
                 cells_to_plot=cells_to_plot,
                 n_genes=10, 
                 colour=example_col,
                 limits = c(-3, 3))

if(DEBUG)cat(file=stderr(), "output$crHeat_plot1:done\n")
p

## $tree_row
## [1] NA
## 
## $tree_col
## [1] NA
## 
## $kmeans
## [1] NA
## 
## $gtable
## TableGrob (5 x 6) "layout": 6 grobs
##   z     cells              name                    grob
## 1 1 (4-4,3-3)            matrix gTree[GRID.gTree.16274]
## 2 2 (4-4,4-4)         row_names   text[GRID.text.16275]
## 3 3 (3-3,3-3)    col_annotation   rect[GRID.rect.16276]
## 4 4 (4-4,2-2)    row_annotation   rect[GRID.rect.16277]
## 5 5 (3-5,6-6) annotation_legend gTree[GRID.gTree.16285]
## 6 6 (3-5,5-5)            legend gTree[GRID.gTree.16288]

if(DEBUG)cat(file=stderr(), "output$crPrioGenes\n")

dt = prioritized_genes[[input$cluster5]]
DT::datatable(dt)

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

print(ls())

##  [1] "Col"                     "DEBUG"                  
##  [3] "allexprs"                "anno"                   
##  [5] "annotation"              "bin"                    
##  [7] "cells.1"                 "cells.2"                
##  [9] "cells_to_plot"           "clustering"             
## [11] "counts"                  "data.1"                 
## [13] "data.2"                  "dataTables"             
## [15] "df"                      "dfOut"                  
## [17] "dge"                     "dt"                     
## [19] "ercc"                    "example_Cols"           
## [21] "example_K"               "example_col"            
## [23] "expression"              "fd"                     
## [25] "featureData"             "featureDataReact"       
## [27] "featuredata"             "filter_by_expr_features"
## [29] "gbm"                     "gbmMat"                 
## [31] "gbm_log"                 "geneC"                  
## [33] "geneNames"               "geneid"                 
## [35] "genes.diff"              "genes.use"              
## [37] "genesin"                 "goodCols"               
## [39] "h"                       "i"                      
## [41] "ids"                     "input"                  
## [43] "ipIDs"                   "kNr"                    
## [45] "km"                      "kmClustering"           
## [47] "line1"                   "line2"                  
## [49] "line3"                   "line5"                  
## [51] "log2cpm"                 "mColor"                 
## [53] "map"                     "maxG"                   
## [55] "medianENSG"              "medianUMI"              
## [57] "merge"                   "mergeExprs"             
## [59] "minG"                    "minGene"                
## [61] "mt"                      "p"                      
## [63] "p1"                      "params"                 
## [65] "pca"                     "pheno_data"             
## [67] "plotexprs"               "positiveCells"          
## [69] "positiveCellsAll"        "prioritized_genes"      
## [71] "rbPal"                   "reads"                  
## [73] "rmIds"                   "rnames"                 
## [75] "scaterReads"             "seed"                   
## [77] "selGenes"                "selectedDge"            
## [79] "selectedgeneList"        "subsetData"             
## [81] "subsetExpression"        "toReturn"               
## [83] "top.genes"               "total.diff"             
## [85] "tsne"                    "tsne.data"              
## [87] "tsneData"                "umiC"                   
## [89] "useCells"                "useGenes"               
## [91] "x"

sessionInfo()

## R version 3.4.3 (2017-11-30)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS Sierra 10.12.6
## 
## Matrix products: default
## BLAS: /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] C
## 
## attached base packages:
## [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] bindrcpp_0.2               scater_1.6.3              
##  [3] SingleCellExperiment_1.0.0 SummarizedExperiment_1.8.1
##  [5] DelayedArray_0.4.1         matrixStats_0.53.1        
##  [7] GenomicRanges_1.30.3       GenomeInfoDb_1.14.0       
##  [9] IRanges_2.12.0             S4Vectors_0.16.0          
## [11] shinyWidgets_0.4.1         kableExtra_0.7.0          
## [13] knitr_1.20                 SCORPIUS_1.0              
## [15] cellrangerRkit_2.0.0       Rmisc_1.5                 
## [17] plyr_1.8.4                 lattice_0.20-35           
## [19] bit64_0.9-7                bit_1.1-12                
## [21] Biobase_2.38.0             BiocGenerics_0.24.0       
## [23] Matrix_1.2-12              reshape_0.8.7             
## [25] mclust_5.4                 RColorBrewer_1.1-2        
## [27] sm_2.2-5.4                 shinyTree_0.2.2           
## [29] threejs_0.3.1              igraph_1.1.2              
## [31] pheatmap_1.0.8             edgeR_3.20.8              
## [33] limma_3.34.9               DT_0.4                    
## [35] shinythemes_1.1.1          plotly_4.7.1              
## [37] ggplot2_2.2.1              shinydashboard_0.6.1      
## [39] shiny_1.0.5               
## 
## loaded via a namespace (and not attached):
##  [1] Rtsne_0.13             ggbeeswarm_0.6.0       colorspace_1.3-2      
##  [4] rjson_0.2.15           rprojroot_1.3-2        XVector_0.18.0        
##  [7] base64enc_0.1-3        AnnotationDbi_1.40.0   xml2_1.2.0            
## [10] ranger_0.9.0           codetools_0.2-15       splines_3.4.3         
## [13] tximport_1.6.0         jsonlite_1.5           readr_1.1.1           
## [16] compiler_3.4.3         httr_1.3.1             backports_1.1.2       
## [19] assertthat_0.2.0       lazyeval_0.2.1         htmltools_0.3.6       
## [22] prettyunits_1.0.2      tools_3.4.3            gtable_0.2.0          
## [25] glue_1.2.0             GenomeInfoDbData_1.0.0 reshape2_1.4.3        
## [28] dplyr_0.7.4            Rcpp_0.12.15           iterators_1.0.9       
## [31] crosstalk_1.0.0        stringr_1.3.0          rvest_0.3.2           
## [34] mime_0.5               irlba_2.3.2            XML_3.98-1.10         
## [37] princurve_1.1-12       zlibbioc_1.24.0        MASS_7.3-49           
## [40] scales_0.5.0           TSP_1.1-5              hms_0.4.1             
## [43] rhdf5_2.22.0           yaml_2.1.16            memoise_1.1.0         
## [46] pbapply_1.3-4          gridExtra_2.3          biomaRt_2.34.2        
## [49] stringi_1.1.6          RSQLite_2.0            foreach_1.4.4         
## [52] rlang_0.2.0            pkgconfig_2.0.1        bitops_1.0-6          
## [55] evaluate_0.10.1        purrr_0.2.4            bindr_0.1             
## [58] labeling_0.3           htmlwidgets_1.0        magrittr_1.5          
## [61] R6_2.2.2               DBI_0.7                pillar_1.2.0          
## [64] fitdistrplus_1.0-9     survival_2.41-3        RCurl_1.95-4.10       
## [67] tibble_1.4.2           rmarkdown_1.8          viridis_0.5.0         
## [70] progress_1.1.2         locfit_1.5-9.1         grid_3.4.3            
## [73] data.table_1.10.4-3    blob_1.1.0             digest_0.6.15         
## [76] xtable_1.8-2           tidyr_0.8.0            httpuv_1.3.5          
## [79] munsell_0.4.3          beeswarm_0.2.3         viridisLite_0.3.0     
## [82] vipor_0.4.5

10X selector report

Bernd Jagla

3/8/2018

parameters